#include "includesT.h"
#include "pocsT.h"
// Code performing Hilbert projection.
using namespace std;
void P1hilbert(ScanData & s, int j, int k)
{
	int i,n, N=0;
	Float pterm=0,x=0, temp =0, temp2=0,a = 0.5 ,b = 1.-a,h;
	vector<Float> image (s.Ny, 0);
	
	
	for (i=0; i<s.Nx;i++){
		if(s.CSmask[i][j][k]) {
			pterm += s.dx*s.recon[i][j][k];
		}
	}
	//cout<< pterm<<endl;
	// Compute H transform


	for(i=0; i<s.Nx;++i){ // x loop
		image[i]=0;
		for(n=0; n<s.Nx;++n){ // x prime loop
			if( abs(i-n)%2 != 0 ) image[i] += s.recon[n][j][k] * (2./(n-i))/PI;
		}
		//s.HT[i][j] = image[i];
	}
	// perform smooting
	temp = image[0];
	for(i=1; i<s.Nx-1;++i){ // x loop
		temp2=image[i];
		image[i] = 0.5*b*temp + a*image[i] + 0.5*b*image[i+1];
		temp=temp2;
	}

	// Compute projection 
	for(i=0; i<s.Nx;++i){
		if (s.DBPmask[i][j][k]) image[i] = s.DBPdata[i][j][k]*(0.5/PI);
	}
	
	// Compute inverse H transform
	for(i=0; i<s.Nx;++i){ // preweight
		x = s.xmin + i*s.dx;
		image[i] = image[i] * sqrt(s.X2 - x*x);
	}
	
	for(i=0; i<s.Nx;++i){ // convolution
		s.recon[i][j][k]=0;
		for(n=0; n<s.Nx;++n){ // x prime loop
			if( abs(i-n)%2 != 0) s.recon[i][j][k] += image[n] * (-2./(n-i));
		}
		
	}
	// Compute preterm


	/* POCS way
	for (i=0; i<s.Nx;i++){
		if(s.CSmask[i][j][k]) pterm += s.dx*s.recon[i][j][k];
	}*/

	// perform smooting
	temp = s.recon[0][j][k];
	for(i=1; i<s.Nx-1;++i){ // x loop
		temp2=s.recon[i][j][k];
		s.recon[i][j][k] = 0.5*b*temp + a*s.recon[i][j][k] + 0.5*b*s.recon[i+1][j][k];
		temp=temp2;
	}
	/* Pterm with outside region
	for(i=0; i<s.Nx;++i){ // x loop
		if(!s.CSmask[i][j][k]){
			N++;
			pterm += s.recon[i][j][k];
		}
	}
	pterm = -pterm/N;
	
	pterm=0;
	for(i=0; i<s.Nx;++i){ // x loop
		
		if(s.APmask[i][j][k]){
			x = s.xmin + i*s.dx;
			N++;
			pterm += PI*sqrt(s.X2 - x*x)*s.APdata[i][j][k] - s.recon[i][j][k];
		}
	}
	pterm = pterm/N;
*/
	for(i=0; i<s.Nx;++i){ // postweight
		x = s.xmin + i*s.dx;
		s.recon[i][j][k] = (pterm + s.recon[i][j][k]) /sqrt(s.X2 - x*x)/PI;
		image[i]=0;
	}
	
}

